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The Parker Scenario for Coronal Heating as an MHD 
Turbulence Problem 
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Abstract. The Parker or field line tangling model of coronal heating is in- 
vestigated through long-time high-resolution simulations of the dynamics of a 
coronal loop in cartesian geometry within the framework of reduced magnetohy- 
drodynamics (RMHD). Slow photospheric motions induce a Poynting flux which 
saturates by driving an anisotropic turbulent cascade dominated by magnetic 
energy and characterized by current sheets elongated along the axial magnetic 
held. Increasing the value of the axial magnetic field different regimes of MHD 
turbulence develop with a bearing on coronal heating rates. In physical space 
magnetic field lines at the scale of convection cells appear only slightly bended in 
agreement with observations of large loops of current (E)UV and X-ray imagers. 



1. Introduction 



Coronal heating is one of the outstanding problems in solar physics. Although 
the correlation of coronal activity with the intensity of photospheric magnetic 
fields seems beyond doubt, and there is large agreement that photospheric 
motions are the source of the energy flux that sustains an active region (~ 
10 7 erg cm~ 2 s~ l ), the debate currently focuses on the physical mechanisms 
responsible for the transport, storage and dissipation (i.e. conversion to heat 
and/or particle acceleration) of this energy from th e photosphe re to the corona. 

A promising model is that proposed by iParkerl ( 19721 . [l988). who suggested 
that coronal heating could be the necessary outcome of an energy flux associated 

with the tangling of coronal field lines by photospheric motions. 

Over the years numerous analytical and numerical investigations ( Parker! 

19721 : I Sturrock and Uchida Ill98ll:lvan Ballegooiienlll986l:lMikic et al.lll989l:|Parkerl 



1988 : lBergerlll99ll:lHevvaerts fc Priest Ul992l:lGomez fc Ferro-Fontan lll992l:lLongcope Sudan 



1994: iHendrix fe Van Hovenl 119961: lEinaudi et all 1996: Georgoul is et alJll99_ 
Dmitruk et al.lll998l : iDmitruk fc Gomez! Il999l : lEinaudi fc Vellil 119991 ) have been 
carried out, discussing in some details different aspects of the problem. 

In all 3D cartesian simulations a complex coronal magnetic field results 
from the photospheric footpoints motions, and though the field does not, strictly 
speaking, evolve through a sequence of static force-free equilibrium states (the 
original Parker hypothesis), magnetic energy nonetheless tends to dominate ki- 
netic energy, and current sheets elongated along the axial direction characterize 
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the system. The results from these studies agreed qualitatively among them- 
selves, in that all simulations display the development of field aligned current 
sheets. However, estimates of the dissipated power and its scaling characteristics 
differed largely, depending on the way in which extrapolations from low to large 
values of the plasma conductivity of the properties such as inertial range power 
law indices were carried out. Furthermore, the physical mechanism responsible 
for current sheets formation was not conclusively determined. 

More recently a first attempt to simulate full 3D sections of the solar coron a 
with a realistic geometry has been performed bv lGudiksen Nordlun d (2005). 
At the moment the very low resolution attainable with this kind of simulations 
does not allow the development of turbulence. The transfer of energy from the 
scale of convection cells ~ 1000 km toward smaller scales is in fact inhibited, 
because the smaller scales are not resolved (their linear resolution is in fact 
~ 500 km). 

While in the future these global simulations will be able to reach the nec- 
essary high resolutions, to investigate the nonlinear dynamics of the Parker 
scenario at relatively high Reynolds numbers, we have recently perfor med high- 



resol ution long-time simulation of the aforementioned cartesian model (jRappazzo et al 
(|2007h ). 

In the next sections we describe the coronal loop model, the simulations we 
have carried out, and give simple scaling arguments to understand the energy 
spectral slopes. 



2. Physical model 

A coronal loop is a closed magnetic structure threaded by a strong axial field, 
with the footpoints rooted in the photosphere. This makes it a strongly anisotropic 
system, as measured by the relative magnitude of the Alfven velocity associated 
with the axial magnetic field va ~ 2000 kms -1 compared to the typical photo- 
spheric velocity u p h ~ 1 kms^ 1 . This means that the relative amplitude of the 
Alfven waves that are launc hed into the corona is very small and, as an efficient 
energy cascade takes place ( Rappazzo et al.ll2007l ). the relative amplitude of the 



fields which develop in the orthogonal planes remains small compared to the 
dominant axial magnetic field. 

We study the loop dynamics in a simplified cartesian geometry, neglecting 
any curvature effect, as a "straightened out" box, with an orthogonal square 
cross section of size t (along which the x-y directions lie), and an axial length 
L (along the z direction) embedded in an axial homogeneous uniform magnetic 
field Bq = Bq e z . This simplified geometry allows us to perform simulations 
with both high numerical resolution and long-time duration. 

The dynamics of a plasma embedded in a st rong axial magnetic field are 
well described by the equations of reduced MHD ( Kadomtsev <fe Pogutsd [l974l : 



Strauss! Il976l ; lMontgome"r"vlll982l ). In this limit the velocity and magnetic fields 



have only perpendicular components, linked to the velocity and magnetic po- 
tentials (p and "0 by 



u± = Vx((pe z ), £>± = V x (V>e 2 ) 



(1) 
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Although numerically we advance the equations for the potentials in order to 
analyze the linear and nonlinear properties of the system it is convenient to write 
the equivalent equations using the Elsasser variables z^ = u± ± b±. The more 
symmetric equations, which explicit the underlying physical processes at work, 
are given in dimensionless form by: 

^ = -(*-. V± )^ + ^ + =W-V ± P (2) 

ot v ' Uph oz 7c 



Vi. • z ± = (4) 

where P = p + b\/2 is the total pressure, and is linked to the nonlinear terms 
by incompressibility @: 



v1p = -£ ^ fe+ . (5) 



hr- 



The gradient operator has only components in the x-y plane perpendicular 
to the axial direction z, and the dynamics in the orthogonal planes is coupled 
to the axial direction through the linear terms tx d z . 

We use a computational box with an aspect ratio of 10, which then spans 

0<x,y<l, 0<z<10. (6) 

The linear terms oc d z are multiplied by the dimensionless parameter v^/uph, 
the ratio between the Alfven velocity associated with the axial magnetic field 
va = Bo/y/4-Kpo, and the photosperic velocity u p h. 

Boundary conditions for our numerical simulations are specified imposing 
the velocity potential <p(x, y) in the bottom (z = 0) and top (z = L) planes: 

y) = i = -. ^2 — ^=p sin I 27r ( kx + 1 V) + 27T £ki] ■ (7) 



These result from the linear combination of large-scale eddies with random am- 
plitudes otki and phases (whose values are included between and 1). We 
excite all the twelve independent modes whose wave-numbers are included in the 

range 3 < (k 2 + Z 2 ) 1 ^ 2 < 4, and then normalize the result so that the velocity 
rms is ~ 1 km s^ 1 . 

In terms of the Elsasser variables z^, to impose a velocity pattern (ttj_) at 
the boundary surfaces means to impose the constraint z + + z~ = 2u*^, and as 
in terms of characteristics (which in this case are simply z^ themselves) we can 
specify only the incoming wave (while the outgoing wave is determined by the 
dynamics inside the computational box), at the top (z = L) and bottom (z = 0) 
planes the following "reflection" takes place: 

z~ = -z + + 2 u° ± at z = (8) 
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z + = -z~ + 2 u\ at z = L (9) 

where and are the forcing functions in the respective boundary surfaces. 

At time t = no perturbation is imposed inside the computational box, i.e. 
b± = u± = 0, and only the axial magnetic field B$ is present: the subsequent 
dynamics are then the effect of the photospheric forcing on the system. 

The linear terms (oc d z ) in equations (f2j)- ([3jl give rise to two distinct wave 
equations for the 2 ± fields, which describe Alfven waves propagating along the 
axial direction z. This wave propagation, which is present during both the lin- 
ear and nonlinear stages, is responsible for the transport of energy at the large 
perpendicular scales from the boundaries (photosphere) into the loop. The non- 
linear terms (z^ ■ Vj_) z^ are then responsible for the transport of this energy 
from the large scales toward the small scales, where energy is finally dissipated, 
i.e. converted to heat and/or particle acceleration. 

An important feature of the nonlinear terms in equations ©-([J]) is the 
absence of self-coupling, i.e. they only couple counterpropagating waves, and if 
one of the two fields z ± were zero, there would be no nonlinear dynamics at all 



This is at the basis of the so-called Alfven effect ( Iroshnikovl 1 19641 ; iKraichnan 



1965), that ultimately renders the nonlinear timescales longer and slows down 



the dynamics. 

From this analysis it is clear that three different timescales are present: ta, 
T p h and T n \. ta = L/va is the crossing time of the Alfven waves along the axial 
direction z, i.e. the time it takes for an Alfven wave to cover the loop length L. 
T ph ~ 5 m is the characteristic time associated with photospheric motions, while 
T n i is the nonlinear timescale. 

For a typical coronal loop ta <C T p h, and for this reason we consider a 
forcing which is constant in time, i.e. for which formally T p h = oo. 

In the RMHD ordering the nonlinear timescale T n i is bigger than the Alfven 
crossing time ta- This ordering is confirmed and maintained troughout our 
numerical simulations. 

The length of a coronal section is taken as the unitary length, but as we 
excite all the wavenumbers between 3 and 4, and the typical convection cell scale 
is ~ 1000 km, this implies that each side of our section is roughly 4000 km long. 
Our grid for the cross-sections has 512x512 grid points, corresponding to ~ 128 2 
points per convective cell, and hence a linear resolution of ~ 8 km. 

Between the top and bottom plate a uniform magnetic field B = Bq e z is 
present. The subsequent evolution is due to the shuffling of the footpoints of 
the magnetic field lines by the photospheric forcing. 

In this section we present the results of a simulation performed with a 
numerical grid with 512x512x200 points, Reynolds number 1Z = 800, and the 
Alfven velocity va = 200 /cms -1 corresponding to a ratio VA/uph = 200. The 
total duration is roughly 500 axial Alfven crossing times (ta = L/va)- 

Plots of the total magnetic and kinetic energies 

E M = ^[dVb 2 ± , E K = ^fdVu 2 ± , (10) 



and of the total magnetic (J) and kinetic (Q) dissipation rates 

irfdvj 2 n = i 
n. J n 



j = L[dVj 2 n = ^-ldvu 2 (ii) 
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Figure 1.: High-resolution simulation with VA/u p h = 200, 512x512x200 grid 
points and 1Z = 800. Left: Magnetic (Em) and kinetic (Ek) energies as a 
function of time (ta, = L/va is the axial Alfvenic crossing time). Right: The 
integrated Poynting flux S dynamically balances the total dissipation D. Inset 
shows a magnification of total dissipation and S for 180 < I/ta < 280. 



along with the incoming energy rate (Poynting flux) S, are shown in Figured) 
At the beginning the system has a linear behavior, characterized by a time linear 
growth rate for the magnetic energy, the Poynting flux and the electric current, 
until time t ~ 6ta, when nonlinearity sets in. The magnetic energy is bigger 
than the kinetic energy, this is the natural result of the field line bending due to 
the photospheric motions both in the linear and nonlinear stages. More formally 
this is a consequence of the fact that, while on the perpendicular magnetic field 
no boundary condition is imposed, the velocity field must approach the imposed 
boundary values at the photosphere both during the linear and nonlinear stages. 

After this time, in the fully nonlinear stage, a statistically steady state is 
reached, in which the Poynting flux, i.e. the energy that is entering the system 
for unitary time, balances on time average the total dissipation rate ( J + O). 
As a result there is no average accumulation of energy in the box, beyond what 
has been accumulated during the linear stage, and a detailed examination of the 
dissipation time series (see inset in Figure [T]) shows that the Poynting flux and 
total dissipations are decorrelated around dissipation peaks. 



Correspondingly, in Fourier space, we have shown (see iRappazzo et al 



20071 ) that the spectral index for total energy fits well the —2 value. We have 



also shown that this spectral index strongly depends on the ratio VA/u p h, i.e. 
on the relative strength of the axial magnetic field. At lower values correspond 
flatter spectra, with an index close to —5/3 (Kolmogorov) , while to higher values 
of the magnetic field the spectra steepens up to ~ —5/2 for VA/u P h ~ 1000. 

Figure [2] shows a snapshot, at time t ~ 18 ta, of current sheets and field 
lines of total magnetic field. 



3. Conclusion and discussion 



Recently a lot of progress has been made in the understanding of MHD tur- 
bulence for a system embedded in a strong magnetic field, both in the con- 
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Figure 2.: Left: Snapshot of current sheets elongated along the axial direction. 
Right: Magnetic field lines are not much entangled because an efficient energy 
cascade takes place. Mid-plane shows field-lines of the orthogonal magnetic field. 
The box has been rescaled for an improved viewing. The box has an aspect ratio 
of 10, to restore the original aspect ratio the box should be stretched 10 times 
along the axial direction z. 



dition of so-called weak (|Ng k Bhattachariee 19971; Goldreich fc Sridhar [l997; 



Galtier et aDteOOOT) and strong t urbulence (iGoldreich k Sridharlll995l.ll9971: | Cho fc Vishniac 



2000; Biskamp k Mulled l2000l: iMiiller et al.ll2003l: iMiiller k Grappinl l2005Tl . In 



particular it has been developed (|Boldvrevi 120051 . 120061 ; iMason et al.l 120061 ) a 
cascade model that self-consistently accounts for current sheets formation. 

The presence of dynamical current sheets elongated along the direction of 
the strong axial magnetic field, that are continuously formed and dissipated, is 
most properly accounted for in the framework of MHD turbulence. In fact, two 
important characteristic of MHD turbulence are that the casca de takes place 
main ly in the plane orthogonal to the local mean magnetic field (jShebalin et al.l 
1983), where small scales form, and that the small scales are not uniformly 
distributed in this plane. Rather they are organized in current- vortex sheets 
aligned along the direction of the lo cal main field, and c onstit u te the diss i pativ e 
structures of MHD turbulence (e.g. iBiskamp k Miillerl (|2000l ). iBiskampI (|2003h 
and references therein). When the system is threaded by a strong axial magnetic 
field, as in our case, the cascade takes place mainly in the orthogonal planes, 
and the current sheets are elongated along the axial direction (Figure [2|). 

The fact that at the large orthogonal scales the Alfven crossing time ta is 
the fastest timescale, and in particular it is smaller than the nonlinear timescale 
r n i (which can be identified with the energy transfer time at the driving scale), 
implies that the Alfven waves that continuously propagate and reflect from the 
boundaries toward the interior are basically equivalent to an anisotropic mag- 
netic forcing function that stirs the fluid, whose orthogonal length is that of the 
convective cells (~ 1000 km) and whose axial length is given by the loop length 
L. 
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The spectra that we have found can be easily derived by order of magnitude 
considerations. Dimensionally, and integrating over the whole box, the energy 
cascade rate may be written as 

e~i\Lp 5 -^-, (12) 

where 5z\ is the rms value of the Elsasser fields z ± = u± ± b± at the perpen- 
dicular scale A. Given the simmetry of the system it is expected and confirmed 
numerically that cross helicity is zero, hence 5z^ ~ 5z^ ~ 5z\. p is the average 
density and T\ is the energy transfer time at the scale A, whic h is greater than 



the eddy turnove r time t\ ~ X/Sz\ because of the Alfven effect (|Iroshnikovll 19641 ; 
lKraichnan1ll965l ). 



In the classical IK case T\ ~ t\ (t\/ta)- More generally, however, as the 
Alfven speed is increased nonlinear interactions become weaker. Simply from 
dimensional considerations as the ratio t\/ta is dimensionless and smaller than 
1, we can suppose that the energy transfer time scales as 

(\ a— 1 
-^J , with a>l, (13) 

where a is the scaling index (note that a = 1 corresponds to standard hydrody- 
namic turbulence). 

Substituting (|13p in (|12p we can compute the energy transfer rate. As this 
is supposed to be constant along the inertial range, considering the injection 
scale A ~ £±, the energy transfer rate (fT2|) becomes 

^ilL.p^^P^ls^ (14) 

To P a n a ~ l 

*-L l ± V A 

On the other hand the energy injection rate is given by the Poynting flux inte- 
grated across the photospheric boundaries: e m = pvA J daup^ • b±. Considering 
that this integral is dominated by energy at the large scales, due to the charac- 
teristics of the forcing function, we can approximate it with 

em ~ p£ 2 ± VAU ph 5z e± , (15) 

where the large scale component of the magnetic field can be replaced with 5zi ± 
because the system is magnetically dominated. 

The last two equations show that the system is self-organized because both 
e and e m depend on 5z£j_, the rms values of the fields at the scale £±: the 
internal dynamics depends on the injection of energy and the injection of energy 
itself depends on the internal dynamics via the boundary forcing. 

In a stationary cascade the injection rate (|15p is equal to the transport rate 
(Q3J). Equating the two determines 6z£ ± , that substituted in (fT4"j) or (fT5"j) yields 
to the energy flux 

02 „„. 2 ( £ -L V A 



£^pv A u p \ h [^) , (16) 



8 



Rappazzo et ah 



where va = -Bo/V^P- This is also the dissipation rate, and h ence the coronal 



heati ng scaling. A dimensional analysis of eqs. ©-(U]) reveals (see lRappazzo et al 
120071 ) that the only free parameter is / = £±va/ Lu p h, so that the scaling in- 
dex a (fl~3|) . upon which the strength of the stationary turbulent regime depends, 
must be a function of / itself, and we have determined its value computationally 
([Rappazzo et al.1 120071 ). 

Dividing eq. (|16p by the surface v\ we obtain the energy flux per unit area 



F = e* Taking for example a coronal loop 40,000 km long, with a number 
density of 10 10 cm -3 , va = 2,000 kms" 1 and u p h = 1 kms , which models 
an active region loop, we obtain F ~ 5 • 10 6 erg cm -2 s" 1 . On the other hand, 
for a coronal loop typical of a quiet Sun region, with a length of 100, 000 km, a 
number density of 10 10 cm" 3 , va = 500 kms" 1 and u p h = 1 kms -1 , we obtain 
F ~ 7 • 10 4 erg cm -2 s" 1 . 

In summary, the Parker scenario for coronal heating may be considered as a 
MHD turbulence problem, self-consistently accounting for current sheets forma- 
tion. The derived coronal heating rates scale with coronal loop and photospheric 
driving parameters, and there is not a single universal scaling law. 
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